Force distribution in a randomly perturbed lattice of identical particles with pair 

interaction 



Andrea Gabrielli 

Istituto dei Sistemi Complessi - CNR, Via dei Taurini 19, 00185-Rome, Italy, 
& SMC-INFM, Physics Department, University "La Sapienza" of Rome, 00185-Rome, Raly 

Thierry Baertschiger 

Dipartimento di Fisica, Universitd "La Sapienza" , 00185-Rome, Raly, 
& Istituto dei Sistemi Complessi - CNR, Via dei Taurini 19, 00185-Rome, Italy 

Michael Joyce 

Laboratoire de Phyisique Nucleaire et Hautes Energies, 
Universite Pierre et Mane Curie - Paris 6, UMR 7585, Paris, F-75005, France 

Bruno Marcos 

Dipartimento di Fisica, Universitd "La Sapienza" , 00185-Rome, Italy, 
& Istituto dei Sistemi Complessi - CNR, Via dei Taurini 19, 00185-Rome, Italy 

Francesco Sylos Labini 

E. Fermi" Center, Via Panisperna 89 A, Compendio del Viminale, 00184 - Rome, Italy 
& Istituto dei Sistemi Complessi - CNR, Via dei Taurini 19, 00185-Rome, Italy 

Abstract 

We study the statistics of the force felt by a particle in the class of spatially correlated distribution 
of identical point-like particles, interacting via a pair force (i.e. gravitational or Coulomb), 

and obtained by randomly perturbing an infinite perfect lattice. In the first part we specify the 
conditions under which the force on a particle is a well defined stochastic quantity. We then study 
the small displacements approximation, giving both the limitations of its validity, and, when it is 
valid, an expression for the force variance. In the second part of the paper we extend to this class of 
particle distributions the method introduced by Chandrasekhar to study the force probability density 
function in the homogeneous Poisson particle distribution. In this way we can derive an approximate 
expression for the probability distribution of the force over the full range of perturbations of the 
lattice, i.e., from very small (compared to the lattice spacing) to very large where the Poisson limit is 
recovered. We show in particular the qualitative change in the large-force tail of the force distribution 
between these two limits. Excellent accuracy of our analytic results is found on detailed comparison 
with results from numerical simulations. These results provide basic statistical information about 
the fluctuations of the interactions (i) of the masses in self-gravitating systems like those encountered 
in the context of cosmological N-body simulations, and (ii) of the charges in the ordered phase of 
the One Component Plasma. 



PACS numbers: 02.50.-r,05.40.-a,61.43.-j 
I. INTRODUCTION 

The study of the statistical properties of the force felt 
by a particle in a gas, and exerted by all the other parti- 
cles of the system through pair interactions, can provide 
useful insights into the thermodynamics or dynamics of 
physical systems in many different contexts and applica- 
tions. Some important examples are given by: (i) the dis- 
tribution of the gravitational force in a spatial distribu- 
tion of point-like masses in cosmological and astrophys- 



ical applications 0, Q, (ii) the statistics of molecular 
and dipolar interactions [J in a gas of particles, (iii) the 
theory of interacting defects in condensed matter physics 
0, and (iv) the contact force distribution in granular 
materials [aQ- 

The seminal work in this field is due to Chandrasekhar 
[H and deals mainly with the probability density function 
(PDF) of the gravitational force in a homogeneous spa- 
tial Poisson distribution of identical point-like masses. 
Specifically the PDF of the force is found to be given 
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exactly by the Holtzmark distribution, which is a three- 
dimensional (3D) fat tailed Levy distribution 0- In 
^ 1^ 1^ approximated generalizations, in different 
physical contexts, of this approach can be found for more 
complex point-like particle systems obtained by perturb- 
ing a homogeneous Poisson particle distribution. 

Here we present a study of the probability distribu- 
tion of the total gravitational (or Coulomb) force for a 
specific class of spatial particle distributions (i.e. point 
processes): three-dimensional shuffled lattices, i.e., lat- 
tices in which each particle is randomly displaced from 
its lattice position with a PDF of the displacement p{u). 
This study can be very useful for applications in both 
solid state physics (e^g. in the case of Coulomb or dipo- 
lar pair interaction) and in astrophysics and cosmol- 
ogy. In the latter context, specifically, large "N-body" 
numerical simulations of self-gravitating particles are an 
essential instrument in the study of the problem of struc- 
ture formation in the universe starting from small initial 
mass fluctuations. These simulations are performed usu- 
ally starting from initial conditions which are suitably 
perturbed simple lattices 11, 12J. We will discuss briefly 
in our conclusions the kinds of questions which may be 
addressed in this context using the formal results and 
methods developed in this paper. 

Quite generally Chandrasekhar ,1] showed that the 
gravitational force acting on each particle can be seen as 
due to the superposition of two different contributions: 
the former (a sort of fluid term) is due to the system as a 
whole, and the latter is due to the influence of the imme- 
diate neighborhood of the particle and therefore depends 
on how a spatial mass distribution (e.g. a fluid) is repre- 
sented through a point-like particle system. The former 
is a smoothly varying function of position while the lat- 
ter is subject to relatively rapid fluctuations. We will 
show that this is true also for the case of a shuffled lat- 
tice, even though some important difference with respect 
to the case of the homogeneous Poisson case are present 
when the shuffling is very small due the extreme uni- 
formity of the particle distribution even at small scales. 
When, instead, the shuffling is sufficiently large we re- 
cover approximately the same behavior as in the Poisson 
case. 

The rest of the paper is organized as follows: 

• Sect. ^ The principal one and two point correla- 
tion properties of point processes and in particular 
of a shuffled lattice are briefly reviewed; 

• Sect, mil We give the statistical definition of the 
global gravitational force acting on each particle 
specifying the conditions under which this is a well 
defined stochastic quantity. In particular we dis- 
cuss the problem posed by the infinite volume limit 
and the necessity of introducing a compensating 
uniform negative background mass density for the 
statistical definiteness of the gravitational force in 
this limit. Here we point out also that this study 



and the definitions given are valid in both the cases 
of gravitational and Coulomb pair interaction; 

• Sect. IIVI We study the stochastic total gravita- 
tional force acting on a particle to linear order in 
the random displacements. In this way we identify 
two different contributions to this force: the for- 
mer comes from the displacements of the point-like 
sources keeping the particle on which we calculate 
the force at its original lattice position, and the 
latter comes from the displacement of this particle 
with respect to the rest of the lattice. The first 
contribution is dominated by the particles in the 
neighborhood of the particle on which we calculate 
the force, while the second can be seen as a force 
due to the system as a continuous mass distribu- 
tion; 

• Sect.[V] Here the variance of the gravitational force 
is calculated in the above approximation and its 
meaning is discussed in relation to the form p{u) 
of the displacement PDF. In particular we find the 
fundamental difference of the behavior between the 
two cases in which the particle displacements are 
limited or not to the elementary lattice cell; 

• Sect. I VII In this section, in order to go beyond 
a study of the statistical properties of the gravi- 
tational force limited to the first and second mo- 
ment, we briefly report some previously known re- 
sults about the PDF of the force in two different sit- 
uations: (i) the flrst is the exact solution of Chan- 
drasekhar for the case of a three-dimensional ho- 
mogeneous Poisson particle distribution, and (ii) 
the second concerns the PDF of the total force in 
one-dimensional (ID) shuflled lattices of particles 
interacting via generic power law pair-interactions; 

• Sect. IVIII We use the results exposed in the pre- 
vious section to develop and discuss an approxi- 
mate evaluation d la Chandrasekhar of the PDF 
(i.e. of all the moments) of the total force acting 
on a given particle. This provides much useful in- 
formation about the exact PDF of the force. The 
approximation becomes better and better when the 
typical permitted displacements of the particles ap- 
proach and go beyond the limit of the elementary 
lattice cell; 

• Sect. rVHll In this section we perform the compar- 
ison of the theoretical and analytical results found 
in the preceding sections with the results obtained 
by numerical simulations. The agreement is shown 
to be very good in general despite the fact that 
the three-dimensional problem of the gravitational 
force in a shuffled lattice is not exactly solvable; 

• Sect. IIXI Finally we summarize the main results 
of this work and draw some concluding remarks on 
their utility. 
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II. STATISTICAL PROPERTIES OF A 
SHUFFLED-LATTICE 

Let us introduce some basic notation that will be useful 
in the rest of the paper. Given a spatial distribution 
of particles in a cubic volume V (we indicate with V 
both the space region in which the system is defined and 
its volume; in this paper we are interested in the limit 
V —^ oo) with equal mass m (i.e. a so called point process 
[TollT^I '). the microscopic mass density function is: 



;^(5(x-x,), 



(1) 



where S{x) is the 3D Dirac delta function, is the po- 
sition of the i*'* particle of the system and the sum is 
over all the particles of the system. Clearly the micro- 
scopic number density is simply given by n(x) = p(x)/m. 
Let us suppose that rip > is the average number den- 
sity. Consequently, the average mass density is simply 
po — nQin. The power spectrum (PS) of such a system is 
defined as 



lim 



{\Shik;V)\') 



nlV 



V^^mm, (2) 



where (...) indicates the average over all the realizations 
of the point process, and 



(5n(k;F) 



is the Fourier integral in the volume V of the number 
density contrast [n(x) — no] HJ. Note that P(k) does 
not depend on m. Therefore without loss of generality 
we fix m = 1 for sake of simplicity and use n(x) for both 
the number and the mass density function. 

In general an initial particle distribution can be per- 
turbed by applying a stochastic displacement to each par- 
ticle of the system (see Fig.^. In particular a perturbed 
lattice is built by applying such a displacement to each 
particle initially belonging to a regular lattice (e.g. sim- 
ple cubic). If R is the generic lattice site the density 
function n(x) for a regular lattice reads 



n(x) 



R 



(5(x - R) 



(3) 



where the sum is over all the lattice sites. If, moreover, 
u(R) is the displacement applied to the particle initially 
at R, the function n(x) for such a SL will be 



n(x ) = 'V 5[x - R - u(R) 



R 



(4) 



In Fig. 121 we present a typical configuration of a 3D per- 
turbed cubic lattice projected on one of the lattice planes. 
The perturbed lattice is said to be uncorrelated if the 
displacements applied to different particles are taken in- 
dependent of each other. We call such a system a shuf- 
fled lattice (SL). This means that p^^) ^ 
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FIG. 1: The figure presents a pictorial view of the effect of 
a stochastic displacement field on a spatial particle distribu- 
tion in 2D. The particles move, through the displacements 
(dashed arrows) , from the old positions (empty circles) to the 
new ones (black circles). 



p[u(R)]p[u(R')], where p^^) [^(r)^ ^(r')] jg the joint 
PDF of the displacements applied to two different parti- 
cles respectively initially at R and R', and p(u) the PDF 
of a single displacement. Without loss of generality in the 
following we limit the discussion to the symmetric case 
in which p(u) = p(— u). Clearly the average over all the 
possible realizations of the displacement field coincides 
with the average over all the possible realizations of the 
point process. Therefore, we will use the notation 



(g(ui,...,u;)) 



d^ui...d^uig{ui, ...,ui)p{ui)...p{ui] 



for the average of any function of the / displacements 
applied respectively to I different points. 

By calling p(k) the characteristic function of a single 
stochastic displacement, i.e., its Fourier transforni (FT), 
we can write the PS of the final point process as [I3 



1 



no 



|p(k)p)-f (27r)3|p(k)p^5(k-H), 

(5) 

where H is the generic vector of the reciprocal lattice 
of the real space lattice giving the initial particle 
configuration, and 



7'in(k) = (2^)^ 



H^O 



H) 



is the PS of this initial regular configuration. A Id ex- 
ample of such a PS is given in Fig. |31 in which the exact 
formula ((Sj is compared with the PS evaluated in numer- 
ical simulations of the SL, showing a perfect agreement. 
Note that if the applied displacement field is isotropic, 
p(u) depends only on u = |u| and j5(k) only on fc = |k|. 
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FIG. 2; The figure provides the projection on the x — y plane 
of a 3D shuffled simple cubic lattice of 16^ particles in a cu- 
bic unitary volume. In this case p(u) — Yl^=i fi'^t) where 
f{ui) — 0{A — |Mi|)/(2A) (as in the simulations considered in 
Sect rTTTTi l with 9.A =f/4. n_e_ each 
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FIG. 3: Comparison between the PS S{k) — n^V{k) measured 
in immerical simulations (circles) of a ID SL with the theoret- 
ical prediction (continuous curve) given by Eq. Q in the case 
in which the one displacement PDF p{u) — 8{A — |it|)/(2A) 
with S = A/£ = 1/50. In order to represent appropriately the 
modulated delta functions contribution to the PS S{k), we 
have normalized their amplitude to a value 10^ for the unper- 
turbed lattice. The numerical result is obtained by averaging 
over 10'' realizations of the same SL and the agreement with 
the theoretical prediction is excellent. 



The connected two-point correlation function [Tol IT^ , 
also called covariance function, is defined in general as 



e(x,x') = 



(n(x)n(x')) 



-1, 



In the case of statistical translational invariance it is a 
function only of the vector distance x — x', i.e., ^(x, x') = 
^(x — x') and satisfies the relation 

e(x)=^-M^(k)]. 

However, as the initial particle configuration is a regular 
lattice, for a SL there is not full translational invariance. 
In this case J^~^[7'(k)] can be seen as the average of 
f(xo -I- X, Xq) over Xq in an elementary lattice cell fs^l- 
Note that for any particle distribution 



e(x) = ^-M^(k)] 



5(x) 
no 



-e(x), 



(6) 



where 5(x)/no is the singular "diagonal part" of the co- 
variance function due to the fact that the microscopic 
density has an infinite discontinuity on any particle, and 
f(x) is the smooth "off-diagonal" part giving the real 
two-point correlation between different particles. One 
important property is ^(x) > —1 for all x. This comes 
from the fact that the average density of particles seen 
by any particle of the system cannot be negative. The 
off-diagonal part ^(x) vanishes identically only for the 
homogeneous Poisson point process [13, 13 • 

It is evident from Eq. ||SJ| that the random shuffling 
{u(R)} in general does not erase completely in the PS 
V(k) the presence of the so-called Bragg peaks (i.e. the 
sum of delta functions) for each reciprocal lattice vec- 
tor H ^ 0, of 7-'jjj(k), but only modulates their ampli- 
tude proportionally to |p(H)p, and adds a continuous 
contribution typical of stochastic particle distributions 
(1 — |p(k)p)/no. Around A; = (more precisely in the 
whole first Brillouin zone [13) 7'in(k) = identically 



(i.e. we can say that around k = 'Pin(k) 



Thus, 



from this point of view, regular lattices can be seen as 
the class of the most superhomogeneous particle distri- 
butions 0,^3' of the particle distributions with a 
PS satisfying P(k) ~ fc" with a > at small k which is 
equivalent to the condition J (Px£,{x.) — 0. As is clear 
from Eq. Q, in this region 7^(k) is determined by only 
the behavior of the displacement characteristic function 
p(k). In particular, even though the lattice is strictly 
anisotropic, this implies that if the displacement field is 
statistically isotropic the final SL has isotropic mass fluc- 
tuations at large scales (i.e. for k ^ 0) . By assuming 
p{u) — p{u), it is now simple to show |l5j that for any 
PDF at small k we have 



p{k) = l-Afc"+o(fc"). 



(7) 



where a = 2 and A = u^/O if is finite (where /(u) = 
/ d^up{u)f{u)), while a = /3 — 3 and A > if p{u) 
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decays as u ^ at large u with /3 < 5 (note that /3 > 3 for 
any PDF p{u) in three dimensions), and therefore is 
infinite. 

In this paper we focus our attention on the class of SL, 
where, as written above, the applied displacements are 
statistically uncorrelated. 



III. DEFINITION OF THE GRAVITATIONAL 
FORCE ON A PARTICLE IN A PERTURBED 
LATTICE 

Let us now consider basic questions concerning the def- 
inition of the gravitational force acting on a particle in an 
infinite perturbed lattice. As aforementioned we assume 
that: (i) all particles have mass m — 1, (ii) the aver- 
age number density is uq = l^^, where (. is the lattice 
spacing, (iii) the microscopic number density is given by 
Eq. and (iv) we choose the units so that the grav- 
itational constant G — \. Let us suppose for simplicity 
that the initial position of the particle on which we eval- 
uate the gravitational field is the origin of coordinates 
(i.e. R = 0). The total gravitational force acting on it 
is: 



R + u(R) -u(0) 
|R + u(R)~u(0)|3 



(8) 



where the sum is over all the particles initially at the lat- 
tice sites R ^ 0. The same sum can be simply rewritten 
as 



(fxn'{x.) 



X - u(0) 



>v |x-u(0)|3 
where V is the system volume, and 

n'(x) = n(x) - ^[x- u(0)] 



(9) 



(10) 



with n(x) given by Eq. Q). 

Note that by simply inverting the sign of the pair in- 
teraction, and therefore of the total force from attrac- 
tive to repulsive, and substituting the identical masses 
with identical charges, Eq. ijHJl gives the total repulsive 
Coulomb force Fc in a perturbed Coulomb lattice 19] on 
the point-charge at u(0) and generated by the all other 
point-charges at R + u(R): 



u(0) - R - u(R) 

''^^w- 



u(R) 



u(0)|= 



u(0) 



U(0)|3 
. 

Consequently, all the results in this paper can be directly 
applied to the statistics of the repulsive force in a shuffled 
Coulomb lattice. 

We are interested in the limit — > oo of Eq. © , where 
we mean by this limit that the volume V tends to the 
whole of K^. It is well known in different physical con- 
texts Tf, "2^] that, if no > 0, the infinite volume limit of 
lattice sums such as those in Eq. © or (|ll|l is in fact not 



well defined because these sums are only conditionally 
convergent in the limit V 00, i.e., their result depends 
on the order in which the single terms are summed. In 
many physical applications, however, as in the case of 
the Coulomb lattice JJJ, this sum is regularized auto- 
matically by the presence in the physical system of a 
uniform background charge density with opposite sign 
with respect to that of the point-charges and such that 
there is overall charge neutrality. Once the attractive 
force Fb of the background on the point-charges is taken 
into account, by adding the corresponding term 



Ff,(uo) =no d^x-^ 
Jv |x 



X- uo 



to Eq. (|ll|l . then the total Coulomb force acting on a 
given point-charge is finite and its value is independent 
of the way in which the infinite volume limit is taken . 
To clarify this point, let us consider the following system: 
a density of point-charges n(x) given by Eq. JQ) with 
m — 1 embedded in a uniform background charge density 
—no — — (n(x)). Therefore the local charge density at 
the point x will be 6n{x) — n(x)— riQ. The force acting on 
a probe charge at the point y of the space and generated 
by the total charge in the volume V around it will be: 



F{y;V) 



d X Sn{'x.)- — 

V \y 



(12) 



Let us now assume that the PS of n(x), and therefore 
of (5n(x), has the behavior 7'(k) ~ fc" at small k, with 
a > — 3 (in three dimensions). Under these hypotheses 
(see e.g. [13) it is simple to show that 



d^a 



6n(x) 



yb/3 



with 6 = 3- aifa<l and 5 = 2 if a > 1. From Eq. IpH) 
we then expect in general that F(y) = F(y; V +00) 
is a well defined stochastic quantity (i.e. a stochastic 
quantity whose statistics is well defined and not depend- 
ing on how the compact volume V is sent to infinity) if 
a > —1 [331 ■ Indeed for these values of a fluctuations 
in the density contrast i5n(x) generate in the infinite vol- 
ume limit quadratic fluctuations in the force F which are 
size independent. This is, in particular, the case for the 
SL we consider in this paper, for which we have shown 
that a > for any displacement PDF p{u) (see Eq. {T))). 
Given that the limit U — > 00 of Eq. (|12|l does not depend 
on the way in which the limit is performed, we can choose 
for simplicity to take the volume V symmetric with re- 
spect to the point y where the force is computed. In such 
a volume the contribution to the force F(y; V) from the 
background vanishes by symmetry. Consequently, with 
this choice of the volume V, the force F(y) coincides with 
the limit of the sum H12|l with 5n(x) replaced by n{x.) (i.e. 
summing only over the particles). 

This in particular implies that for the SL particle dis- 
tribution we will study here, in the presence of an oppo- 
sitely charged neutralising background, we can evaluate 
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the well defined global force F acting on the point-charge 
in u(0) simply using Eq. Hill) where the sum is over all 
the charges contained in a sphere S[r; u(0)] of radius r 
centered on u(0) and then taking the limit r — + oo. Note 
that, on the other hand, the limit r oo of Eq. Hll() 
using instead the sphere S[r; 0] centered on the point 
0, where the considered particle was before the displace- 
ment, does not give the full force F on the point u(0). 
To obtain it the background contribution, which now is 
different from zero, must be added. 

In the rest of this paper, for convenience, we will choose 
to take the infinite volume limit in the former way, i.e., 
symmetrically in spheres about the considered particle, 
also in the treatment of the gravitational problem, i.e., for 
Eq. It is in fact is possible to show that the presence 
of an analogous background with a negative mass density 
(now exerting a repulsive force on the massive particles) 
comes out naturally when the motion of a particle is de- 
scribed in comoving coordinates starting from the exact 
equations of general relativity in a quasi-uniform expand- 
ing Universe (see e.g. HJ). In pure Newtonian gravity, 
instead, such a background does not exist and has to be 
introduced artificially to regularize the problem (Jean's 
swindle, see |22|). The results given here for the statistics 
of the force F have thus to be understood as strictly valid 
in the presence of such a compensating background. 

In |22| some of us (TB and FSL) have given a sim- 
ple estimation of the contribution of the first six nearest 
neighbors (NN) of the particle to Eq. JHJ for a simple cu- 
bic lattice. We show now that instead the sum of Eq. © 
on a sphere of radius r around the central particle can be 
seen as the sum of two different contributions: the first 
"asymmetric" one is due to the self-shuffling u(0) of the 
center particle from the initial R = 0, and can be seen as 
induced by the system as a whole. The second "symmet- 
ric" term is due to the shuffling of all the other particles, 
and is dominated by the contribution of the first six NN. 



IV. THE SMALL DISPLACEMENTS BEHAVIOR 
OF THE FORCE 

In this section we will give the approximate expression 
of F obtained through a linearization in the particle dis- 
placements. The statistical meaning and the limitations 
on taking averages of powers of this linearized expres- 
sion will be discussed in the next section. To simplify 
our computation, but without loss of generality of the 
results, we limit the discussion to those lattices with a 
cubic symmetry: i.e., simple cubic, body centered cubic, 
and face centered cubic lattices. 

As shown above, let us rewrite Eq. ||HJ) as 



are interested in the linear contribution in the displace- 
ment field F^^'^ to Eq. (|13|l for small displacements. Note 
that 



dF(x) = n'(x)- 



:dh 



u(0)P 



is the force contribution coming from the volume element 
d^x around x. Therefore we can rewrite Eq. H13() as 

F = lim F[r;u(0)] = lim / dF(x) (14) 



We now write 



/ dF(x) = [ dF(x) + / dF(x 

^S[r;u(0)] Js{r;0) JsS[r:u{0}] 



) (15) 



where 5(r; 0) is the sphere of radius r centered on the 
lattice point R = 0, and the integration over 6S[r; u(0)] 
means the integration over the portion of S[r; u(0)] not 
included in S{r; 0) minus that one over the portion of 
S{r;0) not included in S'fr ; u(0)] Note that these 

portions of sphere have both the same volume, which 
is proportional to |u(0)|. Consequently, we expect that 
this second integral will give a force contribution of order 
u(0). 

Let us start by evaluating the first term in Eq. (fT^ : 



Fs(r)= / dF(x) 
Js{r-a) 

to first order in the displacements. It can be written as 
Eq. (jSJ where the sum is over all the particles contained 
in the sphere S{r; 0) with the exception of the one at 
u(0). We thus denote this kind of sum by X^r^qi i-^-i 



Fs(r) 



^ IR- 



R + u(R) - u(0) 
u(R) -u(0)|3 



(16) 



We are interested in the contribution F^'"' (r) to this quan- 
tity at linear order in the displacement field. Performing 
a Taylor expansion of Eq. (|16|l to first order in the (rel- 
ative) displacements, we can write: 



R<r 



with 



3[u(0) • R]R- u(0) u(R) 
iR = ^ 



3[u(R) • R]R 



i?3 



(17) 



(18) 



lim / n'(x) , , , ,^ 

— 75[nu(0)] ^ 1X-U(0)|3 



(13) 



where the integral is over the sphere S[r; u(0)] defined 
above, and n'(x) is given by Eq. H1U|) . In this section we 



where R = R/i? is the unit vector in the direction of 

the lattice site R. The quantity F^'-'(r) will be a good 
approximation to Fs{r) if all the (relative) displacements 
are sufficiently small (see next section). Eqs. (|17|) and 
(|18|l show that, at first order, the contribution to the 
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force separates into a part due to the displacement of the 
particle on which we calculate the force, and a part due 
to the displacement of the particle originally at R. Let 
us write 



^R 



(19) 



and we will call simply 



= lim Ff(r) 



with 



.„ _ ~u(0) + 3[u(0) ■ R] 
u(R) ~ 3[u(R) • R]R 



fS 

IR — 



i?3 



(20) 
(21) 



It is simple to show, for the class of lattices we consider, 
that 



R<r 

Rt^O 



(22) 



when the sum is taken up to include all the lattice sites 
R 7^ distributed symmetrically with respect to R = 
up to a given distance. The key points to show this are: 
(i) to rewrite the sum in Eq. 1)22(1 as a sum over all the 
permitted values of i? = |R| < r of the sums over all the 
sites with fixed R; (ii) to consider that for the sub-sum 
over all the sites R, with fixed i? = IRI we can write: 



|R|=iJ 



|R|=i?, 



E fR=-4 E (-u(0) + 3[u(0).R]r) ; 

R R 

(iii) and to note that 

\R\=R |R| = fi |R|=_R 

3 J2 [u(o)-R]R= E (r-r)u(o)= E 

from which Eq. H22|) follows. Therefore we can conclude 
that 

F(0(,)^gu(R)-3[u(R).R]R^ (23) 



R#0 



i?3 



Let us consider now the second term in Eq. (|15|l which 
we call Fyi(r): 



FA{r) - / d 

JSS\r-MO)] 



,3 n'(x)(x-u(0)) 
|x-u(0)|3 



(24) 



Since we wish to evaluate the above integral to linear 
order in u(0), and the volume of integration is already of 
order |u(0)|, we can substitute n'(x) with its average (the 
point u(0) is not included in the volume of integration) 
no = and put u(0) = in the integrand. It is then 
straightforward to show that, taking the limit r oo 
and working to first order in u(0), we have 



An 



nou(O) 



(25) 



where we have simply called Fa — limr^oc Fa (f ) and 
F^^ its approximation at linear order in the displace- 
ment. Note that, as discussed in the previous section, 
this quantity can be seen as the force exerted by the neg- 
ative background contained in the sphere Sir; 0) on the 
particle at u(0). 

Let us summarize before proceeding further in the next 
section. We have now seen that, to first order in the 
displacements, we can write 



J 



r(0 _ 



47r ^ ^ 

— nou(O) -I- lim 



Rt^o 



3[u(R) • R]R 



i?3 



(26) 



r 



where the sum is taken in the sphere S{r;0). As ex- 
plained in Sect.lml Eq. ^ gives the well defined infinite 
volume limit of the force, to first order in the displace- 
ments, generated by the system composed both by the 
particles and the uniform negative background density. 
Note that Eq. (|25|l can be interpreted as either the force 
on the particle at u(0) due only to the other particles 



contained in the sphere 5[r;u(0)] (with r oo), with 
no contribution from the background due to symmetry 
reasons, or the sum of the force on the particle at u(0) 
coming from both the other particles contained in the 
sphere S{r; 0) (with r ~^ oo) and the background in the 
same sphere. 
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V. SMALL DISPLACEMENTS VARIANCE OF 
THE FORCE 

We now turn to the problem of the variance of the 
force (F^)- Considering Eq. (|26|l . one might think that, 

if the variance of the displacements ^ then the 
right hand side of Eq. I|2t)|) can be used to evaluate the 
variance of the force. This is not, however, as we will 
see, always true. In fact, we will show below, if p{u) > 
for M > in such a way to permit at least a pair (and 
therefore an infinite number of pairs) of particles of the 
SL to be found arbitrarily close to one another, then {F^') 
diverges due to the divergence of the pair interaction for 
vanishing separation. This clarifies the meaning and the 



validity of this "small displacement" approximation given 
by Eq. (|26|l : in order to use it to calculate (-F^), it is not 
sufficient to have <C It is instead necessary (and 
sufficient) that all the displacements are smaller than half 
the elementary lattice cell size £ (i.e. that the support 
of p(u) be completely contained in the elementary lattice 
cell) so that each particle has a positive minimal separa- 
tion from its first nearest neighbor. 

In this section we will suppose that this is the case, i.e., 
displacements are limited to a region contained in the 
elementary lattice cell. Given this hypothesis, because of 
the mutual statistical independence of the displacements 
applied to different particles, the average quadratic force 
acting on each particle, to order u^, is then simply: 



r(0|2 



x2_ fi<,. ( |u(R)-3[u(R) -RIR 

—no I -"^ + IhTi 

A / r — >r>r> < ^ 



ii6 



(27) 



If, moreover, we assume that not only the displacements 
of different particles, but also the different components of 
the displacement of a single particle are uncorrelated, it 
is simple to show that, for the class of lattices considered, 
we have 



where Un = V and 



u(R) - 3[u(R) • R]R ^ = 2(|u(R)|2) = 2u2. (28) 
With these hypotheses we can therefore write 



47r 



-no 



2E- 



i?6 



(29) 



where the sum is now over all the lattice vectors R 7^ 0. 
Performing the lattice sum in Eq. (|29|l . for a simple cubic 
lattice, we find: 



E = 12 f -3 j (cinn + C2nn + C3nn 
R#o ^ ^ 



16.11^1 



+ ...) 



(30) 



where Ciun is the relative contribution to the sum 
X^R^o W ^^'^ ^'^^ iisarest neighbor lattice sites of 
the origin R = 0, normalized such that cmn — 1 (giving 
C2nn = 1/4 and cann = 4/81). With these approxima- 
tions, we then find for a simple cubic lattice: 



v/(f5) ^ V'(|F(Op) = .v/(|F(;)p) + (|f(|)|2) = anoUo 

(31) 



a^Wie-l+lyl «5.86. 



(32) 



Hence we can draw a first conclusion: while in the case of 
a homogeneous Poisson particle distribution 1] the grav- 
itational force acting on a given particle is dominated by 
the first nearest neighbor, in the present case it is domi- 
nated by two terms: the former is a global term due 
to the displacement with respect to the rest of the system 
of the particle on which we calculate the force , and the 
latter, F5, is mainly due to the set of the first nearest 
neighbors lattice sites, which all lie at "almost" the same 
distance. As shown explicitly below, because of the sym- 
metries of the initial lattice, the net gravitational force 
in a SL is clearly very depressed, for small displacements, 
with respect to the single nearest neighbor contribution 
of the Poisson case with the same average density no. 



VI. USEFUL RESULTS ON THE PROBABILITY 
DISTRIBUTION OF THE FORCE 

In the next section we will generalize to our case the 
method introduced by Chandrasekhar for calculating 
the PDF P(F) of the gravitational force F in a 3D homo- 
geneous Poisson particle distribution IQJ . As a starting 
point we briefly report in this scetion Chandrasekhar's 
results for this latter case. The specific form of P(F) in 
this case is called the Holtzmark distribution, and for this 
reason we will denote it Ph'(F). Subsequently, we give 
a brief presentation of the exact derivation of the P(F) 
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in a ID SL for a general power law pair interaction. Fi- 
nally, we proceed to generalize Chandrasekhar's method 
to the 3D SL by introducing some ad hoc approximations. 
These results, together with those presented in previous 
sections on the small displacement approximation, will 
provide us with a good qualitative comprehension of the 
behavior of P(F) in a 3D SL when the shape of p(u) 
is varied, and in particular when one passes from dis- 
placements limited to an elementary lattice cell to larger 
maximal displacements. 



A. Gravitational force PDF in a homogeneous 3D 
Poisson particle distribution 

In this subsection we give a brief account of the force 
PDF Ph{F) for a three-dimensional homogeneous Pois- 
son particle distribution which, as aforementioned, is 
called the Holtzmark distribution. The exact mathemat- 
ical derivation of this PDF can be found in [l|. One 
considers a homogeneous Poisson particle distribution in 
a cubic volume with average number density no. Note 
that, because of the total absence of correlations in the 
positions of different particles, the PDF of the force/field 
at a point is the same whether the point is occupied or 
not by a particle. 

While it is not possible to write an explicit expression 
for PniF), is is possible for its FT Ah{^ = J^[Ph(F)], 
i.e., for the characteristic function of F 



with 



AH{^) 



-no- 



4(27rg) 



3/2 



15 



(33) 



Note that this is a typical example of a characteristic 
function of a Levy stable distribution [23l | with exponent 
3/2. The fact that Aniq) depends only on q = |q| means 
that Ph{F) depends only on F = |F| (the particle dis- 
tribution is statistically isotropic). Since ^/^(q) is not 
analytic at q — 0, Ph{F) has a power law tail at large F. 
Specificqlly, as ^^(q) ~ [1 - 4no{2Trqf^^ we have 
Pf/(F) ~ F~^/'^. This shows that moments (F") of order 
a > 3/2 diverge [s^. In particular the variance of the 
force diverges. It is simple to see that this is due to the 
small scale divergence of the pair gravitational interac- 
tion together with the fact that particles can be found 
arbitrarily close to one each other. We then expect the 
same large F scaling behavior of -P(F) for a SL when 
the support of p(u) is larger than the elementary lat- 
tice size, as, just as in the Poisson distribution, one can 
then find an infinite number of pairs of particles with 
members arbitrarily close to one each other. We show 
this below both with the exact results in one dimension 
in Sect. IVI Bl and with the approximate approach a la 
Chandrasekhar in Sec IVIII 

We first recall the limiting behaviors of Ph(P) which 
can be deduced directly from Eq. (|33|l : 

Wh{F) ~ -^F^ for F ^ 0+ (34a) 



Wh{F) 



15 



^3/2^_5/2 f^j. p (-34b) 



F, 



2/3 



and where Wh{F) = A-kF'^Ph{"P) is the PDF of the mod- 
ulus of F. The quantity can be seen as the typical 
force acting on a particle and is called the normal field in 
0. It is also important to note that, in the Poisson case 
for large values of F, the Wh{F) is well approximated 
by the PDF of the modulus of the force due only to the 
NN particle [13: 



Wnn{F) ^ 27rnoF-5/2 pxp 



47rrtoF^ 



-3/2 



(35) 



This means that the main contribution to the force felt by 
a particle in a homogeneous Poisson distribution comes 
from the first NN particle, and is due to the small dis- 
tance divergence of the pair gravitational interaction. 



B. Exact results for the ID SL 

Before introducing an approximate approach d la 
Chandrasekhar for the SL in 3D, we give some exact 
results obtained for the force PDF in an analogous ID 
SL of particles interacting via a power law interaction as 
presented in [23 |. 

Let us consider a ID SL, i.e., a set of 2N + 1 point- 
like particle of unitary mass placed at the points Xm = 
m£+Um with m — —N, N of the segment [— L/2, L/2], 
where £ — L/{2N+1) is the lattice spacing (and uq — l/£ 
the average particle density) , and Um is the displacement 
of the TO*'* particle from the lattice position m£. We 
assume that all the Um are extracted from the same PDF 
p{u) independently of one another, and that the particles 
interact via the attractive pair force: 



/(^) 



\a + l 



where x is the particle separation and a > 0. 

Therefore, the particle at xq = uq feels the total force: 



F 



mi + U„i — Uq 



(36) 



Note that F is not a sum of statistically independent 
terms as each of these terms depend on two displace- 
ments, one of which is always uq. By considering that 
the PDF of x„i is simply given by p{xm — rn£), we can 
formally write the PDF P(F) of the stochastic force F 
as i2Jj 



10 



P{F) = 



+ 00 



+ CXD 



Y[ dXmP{Xm~mt) 



— 00, +00 



m/0 



(37) 



r 



in which we have taken directly the (symmetric) Umit one can write the characteristic function of the stochastic 
N ^ 00 with £ fixed. By taking the FT in F of Eq. (|37|), force as 



dxopixo) 



~00, + C30 

n 



dsp{s + Xq — n£) exp 



iqs 



(38) 



r 



Through an exact analysis |2^] of the small q behavior 
of Eq. one can distinguish essentially two classes of 
SL for what concerns the large F behavior of P{F): 

(I) First class: No particle can be found arbitrarily 
close to any other particle; i.e., the supports respectively 
of p{u) and of p{u — n£), for all integer n ^ 0, have 
an empty overlap. Specifically this is the case if 3A e 
{0,£/2) such that p{u) — for > A, i.e., when the 



support of p{u) is all contained in one elementary lattice 
cell. In this case P{F) vanishes at large F for all values 
a > faster than any negative power of F and all the 
moments {F") are finite for any n > 0. If moreover 
p{u) = p{—u), P{F) is not far from a Gaussian with 
zero mean, even though there are deviations from pure 
Gaussianity depending on the exact shape of p{u) . The 
finite variance of F in this case is given by: 



+00 +00 



+00 



1=1 n=l 



+ ({{u -X0 + n£)-'-)^ {{u ~xo + n£)-^)l)^ 



(39) 



r 



where simply ((...))„ = dup{u)[...) and ((...))^^ = 

dxop{xo){...). Thus the force variance is composed 
of two distinct contributions: (i) the double sum which is 
determined basically by the fluctuations created by the 
stochastic displacement xq of the particle initially at the 
origin on which we evaluate the force (in this term the 
average over u, i.e., over the sources, plays the role of a 
smoothing); (ii) the single sum which is instead mainly 
due to the fluctuations in the displacements u of all the 
sources of the force (in this term only the average over 
Xq is a smoothing operation). This coincides qualita- 
tively with what we have seen in Eq. (|29|l for the three 
dimensional case with an approximate calculation. The 
analogy with Eq. (|29|l can be made stronger with the hy- 
pothesis of small displacements, i.e., u'^/£'^ <C 1 where, 
in analogy with 3D, we have called ~ I-^ dup{u)u^. 
Keeping only terms up to the second order in the random 



displacements, it is simple to show that Eq. H39|) can be 
rewritten as: 



-C(2a + 2) 



(40) 



where (^(x) is the Riemann zeta function. In this expres- 
sion the flrst contribution comes from the fluctuations 
of the position of the particle on which we calculate the 
force, and the latter comes from the fluctuations of the 
positions of the sources. In particular by writing 



C{2a 



E: 

n=l 



1 



,2a+2 ' 



it is simple to show that the n^^ term in the sum gives 
the relative contribution of the n*'' NN lattice sites of the 
origin to the force for u^/f^ <C 1. 
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(II) Second class: At least one pair of particles (and 
therefore an infinite number of pairs) can be found arbi- 
trarily close to one each other; i.e., the supports respec- 
tively of p{u) and of p{u — na), for at least an integer 
n ^ 0, have a finite intersection. The simplest case in 
this class is when 3e > such that p{u) > for all 
|u| < £/2 + e. In this case it is possible to show 
that the large F tail of P{F) is proportional to p-^-'^f". 
More precisely 



P{F) ~ BF^i^i/" for F 



(41) 



with 



B = - I dxop{xa) p{xo-ne), (42) 



-2«* /e<n<2u' 



where u* is such that p{u) > for u < u* and vanishes for 
u > u* (and u* — -\-oo for p{u) with unlimited support). 
Note that the large F exponent of P{F) is independent 
of the details of p{u). Moreover it coincides with the 
exponent characterizing the large F tail of the ID Levy 
PDF found for a ID Poisson particle distribution [33, 
i.e., the ID analog of the Holtzmark distribution seen in 
the previous section. In this case the amplitude B of the 
tail in the SL is £ dx^ p{xo) E-tu* /^<r,<2«* ^(^^o - 
n£) < 1 times smaller than the amplitude of the Poisson 
case which is simply l/{a£) 24] . We note also that, if 
u* ^ £ and p{u) is smooth (i.e. approximately constant) 
on the length scale £, we can approximate Eq. (|42|) with 



B 



1 - p{0) £ 
a£ 



(43) 



This last approximated expression is again independent 
of the details of p{u) for u 0. 

Intermediate cases between (I) — rapidly decreas- 
ing P{F) — and (II) — Holtzmark-like power law tail of 
P{F) — are possible only if displacements are permitted 
exactly up to |m| = £/2 but not beyond this value. In 
this case, depending on how p{u) behaves in the neigh- 
borhood oi u = ±£/2, one can have intermediate large 
F tail behaviors of P{F), e.g., a power law decreasing 
faster than i^-i-V". 

The results just outlined coincide qualitatively with 
those we will now find using an approximate generaliza- 
tion of Chandrasekhar's approach to the case of a three- 
dimensional SL. We will see that this method gives accu- 
rate predictions on the large F behavior of P{F) for all 
p{u), even though the accuracy for the amplitude of this 
tail is good only in the limit of sufficiently large displace- 
ments. 



VII. APPROXIMATE CHANDRASEKHAR 
APPROACH TO P(F) IN THE 3D SL 

We extend here the formalism developed in j3| in a sim- 
ilar way to what has been done in 9] for particle distribu- 
tions generated by a Gauss-Poisson process. As shown in 



this latter paper, for spatially correlated point processes 
in which each point has the same mass, it is possible to in- 
troduce an approximated PDF for the gravitational (or 
Coulomb), force felt by each particle (identical charge) 
of the system and due to all the others. The approxi- 
mation consists in using only the information carried by 
one and two-point correlation functions of the number 
density field. 

Let us consider the general problem of a statistically 
homogeneous particle distribution in a cubic volume of 
side L (and hence of volume V = L^) and mean number 

— 1/3 

density tiq (with L ':$> ). As usual the microscopic 

density is given by Eq. with m = 1. We study the 
PDF of the total gravitational force on the particle at 
Xq due to the other N points in the system in the limit 
V ^ oo with N/V = no > fixed (taking, as explained 
above, the limit symmetrically with respect to the point 
xo) llg- For simplicity let us take a coordinate system 
such that Xq coincides with the origin. The other N 
particles occupy the positions Xi,X2, ...,XAr respectively. 
The force acting on the probe particle at the origin is 



N ^ 



(44) 



Let pi^''(xi,X2, ...,xjv) be the conditional joint PDF of 
the positions of the N other particles with respect to the 
probe at the origin. The approximation we use consists 
essentially in assuming the validity of the following fac- 
torization 



(xi,X2, ...,Xjv) = ]Jpc(Xj) , 



(45) 



where Pc(x) is the conditional PDF of the position x 
of a given particle with the condition that the origin of 
coordinates is occupied by another particle of the sys- 
tem. We thus approximate the system seen by the parti- 
cle at the origin with an inhomogeneous Poisson particle 
distribution with space dependent average density pro- 
portional to Pc(x). This hypothesis works well in the 
Gauss-Poisson case jOj, and we expect it not to be bad 
for any particle distribution which does not differ too 
much from a Poisson one, i.e., with a two-point correla- 
tion function ^(x) which is short range (i.e. integrable) 
and with a small amplitude. This means that in our case 
of a SL we expect that the approximation will give a 
quantitatively good estimate of P(F) only when the typ- 
ical displacement of a particle starts to be of the order 
of the lattice cell size. In fact it is only in this case that 
the lattice Bragg peaks contribution to the PS of Eq. (0 
is strongly reduced and, consequently, the amplitude of 
^(x) is small enough. However we will see that even for 
smaller displacements, when the force variance is finite, 
this approximation gives good quantitative predictions 
about the large F tail of ^'(F), even though the value of 
the force variance is accurate only for the case in which 
the largest permitted displacement starts to approach the 
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cell boundary. This means that when instead the max- 
imal permitted displacement is much smaller than the 
lattice cell size we keep the qualitative results we present 
in this section but use for the variance of F Eqs. (|27|) and 
(123. 

Directly from the definition of the two-point correla- 
tion function, we have 

p,(x)=A[l + e(x)], (46) 

where A ~ 1/V is the normalization constant, and ^(x) 



for any stationary point process is defined as the off- 
diagonal part of the covariance function ^(x) [see Eq. ©]. 
In fact (n(x))p — no[l + gives the average condi- 

tional density of particles seen by point occupied by a 
system particle ,1Q.| (where indicates a conditional 
average) . 



Making these hypotheses, it is possible to write [ij 
in the infinite volume limit, the PDF of F as 



P{F)^ J (03e^'' ''exp|-no / d^x [1 + e(x)](l - e-'''-/-' 



It will be useful to rewrite this in the form 

d3 



P(F) 



i2n) 



— e*i-^ exp 



-nQCH{q)-no / (PxCix.){l 



3 — iq-x/a; 



(47) 



(48) 



where Cniq) — 4(27r(7)'^/^/15, with the multiplicative 
factor — ng, is the cumulant generating function for F 
in the homogeneous Poisson case already considered in 
Sect. Wn\ [see Eq. The function 

A(q) = y d3i^e^"i-^F(F) (49) 

= exp I -no y" d^x [1 +e(x)](l - e-*i-^/^')| 

is the characteristic function of the stochastic force F. 
We recall that the function 

5(q) =lnA(q) - -no ^ [1 + e(x)](l - e-"!'"/"') 

. 

is the cumulant generating function of the stochastic field 
F [2^. The cumulants (i.e. the connected parts of the 
moments) of F can be directly calculated by taking the 
derivatives of this function at g = 0. Therefore the small 
q behavior of ^(q) describes the large F properties of 
P(F). Since in Eq. (|49|l the small q region corresponds 



to the small x region of ^(x) we can say roughly that 
the large F behavior of P(F) is basically determined by 
the small separation properties of the particle distribu- 
tion [and therefore on the small x behavior of ^ (x)] . We 
have already considered this aspect both for the Chan- 
drasekhar method for the homogeneous Poisson particle 
distribution and for the exact results in ID. For the 
homogeneous Poisson case the off-diagonal part of the 
reduced correlation function is ^(x) = 0, and Eq. H49() 
becomes consistently Eq. (I33|l which implies Eq. (|34a|l . 
A. The large F behavior of P(F) 

In order to simplify the calculations which follow, we 
make the assumption that ^(x) = ^{x) even though, rig- 
orously speaking, this is not the case for our SL because 
of the underlying lattice symmetry even when p{u) de- 
pends only on u = juj. For the SL this can be seen as an 
approximation consisting in substituting ^(x) in Eq. H48|) 
with its angular average. Assuming ^(x) — ^(x) implies 
that A{q) and P(F) depend respectively only on q and 
F, and Eq. (|50|l can be rewritten as 



g{q) = -47rno / dx x^ [1 + ^{x)] 



-no { Cniq) 

r 



An 1 dx X ^(x) 
'0 



X 

— sm I 

q \x^ 



(51) 



Let us now analyze how the shape of ^ {x) determines the 
large F tail of P(F). To do this the fundamental step is 
to study the small q behavior of Q (q) . In this respect we 



distinguish below three cases for the choice of p(u) with 
a continuous and convex support: (Sect. IVII A II — large 
displacements) it permits, at least in some directions, 
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displacements of particles beyond the border of their ele- 
mentary lattice cells, allowing, in this way, different parti- 
cles to be found arbitrarily close to one each other; (Sect. 
IVII A 21 — marginal displacements) it permits, at least in 
some directions, displacements exactly up to the border 
of the elementary lattice cell and in no direction beyond 
it, allowing different particles to be found arbitrarily close 
to one each other, but only marginally; (Sect. IVII A 31 
— small displacements) its support is all contained in a 
internal region of the elementary lattice cell so that there 
is a finite lower bound on the distance between any two 
particles. 

As has been noted in the discussion above about the 
accuracy of the approximation (|45|l and H46|l . we expect 
to obtain better and better approximations for P{F) the 
larger is the typical particle displacement. 

In all three cases, it is important to note that G{0) — 0. 
Moreover, as a SL is "superhomogeneous" , / d^x£,{x.) — 
—uq — We now treat one by one the above three 

cases. 



1. Large displacements 

In this case it simple to show that ^(0) > —1 and 
finite. In fact the average conditional density no[l + ^(x)] 
has to converge to a positive constant for a: ^ 0, as 
the large displacements permit couples of particles to be 
found arbitrarily close to one another. This is sufficient 
to show (see Appendix that, up to the leading term 
at small q, we have 



^(q) 



15 



no [l + e(0)] (2719)3/2^ 



which is of the same order in q of Ch (q) ■ By recalling that 
A(q) — exp[t/(q)] and performing the Fourier transform 
(I48|l . we can simply derive at large F that 

F(F)~ [l+e(0)]PH(F), 

which can be rewritten in terms of the PDF of F = |F|, 
W{F) = 47rF2p(F), as 



WiF) [l + mWH{F) ~ [i+m] 



15 



p3/2p-5/2 



where we have used Eq. Ij34b|) . This means that in this 
case P{F) presents the same large force scaling behavior 
as that of the Holtzmark distribution, but with an ampli- 
tude greater by a factor 1 + ^(0). Given that no[l -|- ^(0)] 
is the average conditional density of other particles at 
zero distance, it is simple to show that for a SL one can 
write 

1 + ^(0) = — V / d\p{u)p{u ~ R) , 

from which one finds the 3D analog of Eq. H42(l . Note that 
for most choices of p(u) one has also ^(0) < (i.e. the 



system is negatively correlated at small scales). There- 
fore the amplitude of the tail will usually be reduced with 
respect to that in the Holtzmark case. In any case, as for 
the Holtzmark distribution, all the generalized moments 
(F"^) diverge for ^ > 3/2. 



2. Marginal displacements 

In this case ^(x) ~ [-1 -I- Bx^] with B, /3 > at smafi 
X. In fact, as displacements are permitted up to exactly 
the cell boundary, the probability of finding a particle 
at distance x from another fixed particle must vanish 
continuously for x 0. By studying the small scale 
behavior of t/(q) (see Appendix IXji . we can conclude that 
for g — > we have: 

• e(q) - -q(3+/5)/2 

if < /3 < 1, which implies WiF) - F-(5+/5)/2 
at large F. In this case the variance {F^') is thus 
divergent, but slower than in the Holtzmark case; 

• 5(q) ^ q^^ogq 

if /3 = 1, implying W{F) ^ F^^ (giving a logarith- 
mically divergent variance (^F^')); 

• ^(q) - -q^ + o{q^) , 

for j3 > \ where o(g^) is a power vanishing faster 
than 2 and including the main singular part of this 
small q expansion proportional to q(3-i-/3)/2 (■■^^j^j^ 
logarithmic corrections for (3 integer larger than 1). 
This implies again W{F) - F-^^+'^^Z^ g^^- j^^j-gg p 
(giving a finite variance {F"^)). 

Summarizing, we can say that in general W{F) ^ 
p-{5+0)/2 ^j^j^ /3 > at large F implying that aU the 
moments (F'') diverge for ?7 > (3 -I- /3)/2. 



3. Small displacements 

In this last case displacements are permitted up to a 
distance A < £/2. Consequently, there will be a positive 
minimal distance x* = £ — 2A between any two particles. 
This implies that ^(x) = — 1 identically for x < x*, as 
around any particle there is a minimal empty region of 
radius x* . Therefore Eq. 151|) can be rewritten as 



/•OO 

CJ(q) = -47rno / dxx'^[l 

J X* 



1 sm — 

q \x^ 



(53) 

The large F behavior of Eq. 148|l is essentially determined 
by the small q behavior of this function. Since x* > 0, 
and ^(x) for x oo, one can verify easily that the 
integral in Eq. (|53|l can be expanded in Taylor series to 
all orders in q with finite coefficients, i.e., G{(1) and, con- 
sequently, A{q) are analytic functions of q. From Fourier 
transform theorems [2^, one can then infer that P{F) 
vanishes at large F faster than any negative power of F, 
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i.e., W{F) has all moments (F") finite. To the leading 
order in we can write 



1 



■ (^) 



6 



(54) 



Therefore at small q the characteristic function A{q) has 
the following behavior: 



where 



A(q) = 1 



da; 



0(<7*) 



l+e(a^) 



(55) 



(56) 



is the approximation we obtain with this method for the 
variance (-F^) of the force F. It is important to stress 
that this formula for a^p is expected to apply to our SL 
only when IS./ 1 approaches sufficiently the value 1/2, i.e., 
when displacements are large enough to make the ampli- 
tude of ^(x) small. For smaller values of A, Eqs. H27() 
and H29|l have instead to be applied to calculate {F'^)- 

Note that the 3D isotropic (i.e. monovariate) and un- 
correlated Gaussian PDF reads 



Pom 



and has a characteristic function 



3 


3/2 


' 3F2\ 


_27rcr| 


exp y 





^G(q) = exp 



(57) 



(58) 



which has the same second order small q expansion as 
that of Eq. H55|l [i.e. the same lowest order cumulant 
generating function t/(q)]. Therefore we can say that in 
our case, if A is not much smaller than £/2, P(F) is ap- 
proximately given by the 3D Gaussian lf57|) where (i^^) 
is approximated by a^. which is given by the covariance 
function via Eq. (|56|) . Note, however, that Gaussianity is 
only approximate. In fact terms of order equal or higher 
than q^ from the expansion of Eq. (|53|l are in general 
not vanishing differently from the case of Eq. (|57|) where 
the cumulant generating function is a simple quadratic 
function of q. Instead, the fourth order term of the Tay- 
lor expansion of Eq. H53|l can be written in general as 
47rnoC4g^ (the quantity (7r/6)noC4 is the fourth cumulant 
of F) with 



C4 



ax r; 



In order to evaluate how large is the deviation from pure 
Gaussianity due to the fourth (and higher order) cumu- 
lant term, one has to compare AirtiQCi with /72. In fact 
the fourth order term of the Taylor expansion of the pure 



Gaussian Aaiq) (Eq. ^) is {(j%/72)q^, while m our 
case the fourth order term of A{q) is (cr|,/72 + 47rnoC4)g^. 
Therefore the quantity 



Aatg = 



2887rnoC4 _ 3 j^. dx[l + i{x)]/ x^ 
a% ^ 207rno dx[\ + i{x)]/x^' 



(59) 



is a good measure of the degree of the non- Gaussianity of 
P(F). The quantity AjvG is called in probability theory 
the kurtosis excess |27l |. It measures the importance of 
the large F tail with respect to the Gaussian case with the 
same variance. When \ng *C 1 (i-e- 2A/^ ^ 1) we can 
say that the deviation of P(F) from the Gaussian Pg(F) 
in Eq. (|57|l is small, while, on the other hand, if instead 
Xmg — li the deviation starts to be appreciable and the 
large F tail of P(F) starts to be considerably fatter than 
that of Pi/(F), and finally for A^vg ^ 1 (i-e. for A ~ 
£/2) the Gaussian approximation is inappropriate and 
P(F) starts to develop the power law tail described above 
for the cases of large and marginal displacements. This 
deviation from Gaussianity (see below) clearly increases 
with A and in general diverges when A approaches £/2: 
in fact for this value all the moments higher than a given 
value diverge. Therefore for 2A l~ we expect to see 
large discrepancies of P(F) with respect to Eq. (|57|l . 

In general for small (x~ x*^ > the covariance ^(x) in 
the present case behaves as [1 -f i,(x^ = B{x — x*)^ with 
positive B and depending on p(u). By changing the 
value of /3 (i.e. in our case of a SL, by changing the scahng 
behavior of p(u) in the neighborhood of tt = A), we 
can have or not a diverging behavior of a\ and C4 when 
X* — > 0+, that is when A {l/2)~ , and displacements 
are permitted up to exactly the lattice cell boundaries. 
However in general we can say that the deviation from 
pure Gaussianity given by Aatg increases with A and 
diverges at A = £/2 for /3 < 5. 

For what concerns the approximate variance cr|n, it is 
simple to show that if (i) /3 < 1 this quantity diverges 
as x*t^'^ - (£ - 2A)^"i, if (ii) /3 = 1 it diverges as 
— \og{l — 2 A) (as we have already seen above in the case 
in which A = ^/2 exactly), and if (iii) /3 > 1 it converges 
to a finite value. In a similar way it is simple to see that 
C4 diverges as (^ — 2A)^~'^ for /3 < 5 (implying a large 
F tail slower than 1/P^, for P(F)), logarithmically for 
(3 — h, and converges otherwise. 

Finally Gaussianity (i.e. a PDF given by Eq. ^bl\ ) 
is almost exact when A <C £/2 [s^. However in this 
case the approximation given by Eq. I|56|l for the variance 
(P^) of F in the SL is not a good one. In fact in this 
case substituting such a SL with inhomogeneous Poisson 
particle distribution with a radial density around the ori- 
gin (where we calculate the force) equal to the average 
conditional density of the SL is a bad approximation as 
^(x) acquires large values around the Bragg peaks. Nev- 
ertheless, following also the results in d = 1 presented 
above, in this case we can say again that F is approxi- 
mately a 3D Gaussian variable [i.e. with a PDF given by 
Eq. EZI)], but with (P2) given by Eq. 
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B. Small F behavior of P(F) 

In order to find the small F behavior of P{F), first of 
all we note [see the first formula of Eq. H34a|) ] that in the 
homogeneous Poisson case 



is finite, i.e., Wh{F) ^ F"^. In our case, from Eqs. 
and if^ . P(F) for F = is given by 



PniO) - 47r / dgq^e^vi^noCniq)] 



J 



P(0) 



J d^qA{q) = 47r y dq q^ exp 



-noCniq) 

r 



Atttiq j dxx ^{x) 
'o 



X . f q 
— sm 



(60) 



It is simple to verify that for any possible covariance 
functions S,{x) the quantity P(0) stays finite, i.e., again 
W{F) ~ F"^. Roughly speaking, the more the parti- 
cle distribution shows anti-correlations, the larger will 
be the value P(0), i.e., the larger will be the probabil- 
ity of observing a small value of F . On the contrary, 
the larger the positive correlations the smaller the value 
of P(0). In particular, in our case of a randomly per- 
turbed lattice the system is superhomogeneous, that is 
/ d'^x£,{x) = — 1, and therefore negative density-density 
correlations are always more important than positive cor- 
relations. This means that in general, given the structure 
of Eq. I|60() (and in particular the fact that sint/t < 1 for 
any t > 0), P(0) > Ph{0)- Only in the limit of random 
displacements in the whole system volume do we have 
P(0) Pff (0). Moreover in general, the smaller are the 
typical displacements, the larger will be the contribution 
of anti-correlations to Eq. H6U|) and then the larger P(0). 



the case in which the three components Ui with i — 1,2,3 
along the three orthogonal axes of the displacement u ap- 
plied to the generic particle are statistically independent 
and uniformly distributed in a symmetric interval, i.e.. 



(61) 



with 



fiu^) 



4=1 



2^ for \u,\ < A 



(62) 



otherwise . 
In this case the average quadratic displacement is 



- 3 



VIII. COMPARISON WITH NUMERICAL 
SIMULATIONS 



or I P — 6"^ in normalized units 5 ~ Ajl. The FT 
of p(u), i.e., the characteristic function of the random 
displacements, is simply given by 



In this section we compare our theoretical predictions 
for the statistics of F (specifically, the variance and the 
global PDF) we have given in the previous sections with 
numerical results for the same quantities obtained di- 
rectly by computer simulations of the SL particle distri- 
bution with given p(u). The paradigmatic example, on 
which we concentrate our numerical analysis, is given by 



5(k)=.Fb(u)]=n 



sin(A;j A) 



where kj is the j*'' component of k. By using Eq. Q 
it is then simple to verify that the PS 'P(k) of the SL is 
given by: 



J 



P(k) = ^3 



_ A sin^(fcj-A) 
y (%A)2 



A sin^(gj A) 



(5(k-H), 



(63) 



r 



This expression reduces at small k to 0, 0| 



p(k) = \e'A^k^ = \iH''k^ , 



which depends only on k = |k|, i.e., mass (or number) 
fluctuations are statistically isotropic on large scales even 



(64) 
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though the SL is not because of the underlying lattice 
symmetry. In Fig. O the PS of the Id analog of such a 
SL with S < 1/2 is given. In this figure both the fc^ 
scaling behavior at small k and the modulation of the 
Bragg peaks of the initial lattice by the factor are 
clear. 



For our particular choice of p(u) it is possible to calcu- 
late the inverse FT of Eq. (|^ to obtain, through Eq. 
the off-diagonal covariance function ^ (x) for all values of 
A. Let us call n the integer part of the ratio 4A/£. We 
can write 



l+e(x) = 



2A 



3 3 

n 



\ 2A 



1 j 6i(2A- 



^.i)+n 



A^ 

'a) ^(2A)2 ^ 

I 



+ CX2 



-m£\)0{2A' -\xj -me\) , 

m 



where 9{x) is the usual Heaviside step function, and 
2 A' = (2 A - ni/2) if n is even or zero and 2 A' = 
[{n+l)£/2 — 2A] if rt is odd. It is important to notice that 
< 2A' < i/2 for all A values. The function ^(x) given 
by Eq. H65(l has in general a very complicated oscillating 
form, with the exception, as shown below, of the case in 
which A = mi/2 exactly with m integer. However for a 
generic A we can say that it is continuous and is com- 
posed of two different contributions. The former is given 
by the first product of Eq. H65|l and comes from the con- 
tinuous (i.e. purely stochastic) part of the PS, and the 
latter, given by the second product which is a lattice pe- 
riodic function, comes from the modulated Bragg peaks 
part of the PS. For all the choices A = mi/ 2 with inte- 
ger m the "Bragg peaks contribution" exactly vanishes 
leaving only the first stochastic contribution: this means 
that for these values of A the system becomes statistically 
translationally invariant. Note that for A < £/2 we have 
^(x) = — 1 identically in the cube of side \xi\ < {i — 2 A) 
for all i — 1,2,3, meaning that the probability of find- 
ing a particle in such a cube centered around any other 
particle is strictly zero for these values of A. 

We can now draw the following general conclusions for 
all the possible choices of 5 = A/i. For simplicity we 
start with the limiting case 6 = 1/2, and then we analyze 
respectively the cases 6 < 1/2 and S > 1/2. 



A. Large F prediction for A — i/2 

For S — 1/2, it is straightforward to find the exact 
result 

3 

a^) = -l[{l-\x,\/i)e{i-\x,\). (66) 
i=i 

Thus ^(x) is in this case non- vanishing, and negative, 
only in the cube —i < xj < i for all j — 1,2,3. As 
mentioned above in this case the lattice Bragg peaks of 
Eq. H63|) are completely erased by the displacements and 
the system is statistically invariant under translations. 



Expression at x <^ i gives 



?(x) 



3 

E 



(67) 



Therefore the SL falls in the class of Sect. IVII A 21 
with /3 = 1, i.e., with P(F) F'^ at large F [or 
equivalentlyVF(F) ^ F^"^] and, from Eq. ifSBI) . loga- 
rithmically diverging variance ap. In order to have a 
more quantitative description of this case, we mimic the 
anisotropic Eq. (|66|l with the following isotropic ^{x): 



-1 



TTN 1/3 X 

6/ I 



1/3 



i^ 



(68) 



i.e., an isotropic function with a linearly increasing be- 
havior, similar to the one in Eq. (|67|l . from ^' = — 1 at 
x = to = at the border and outside the sphere 
centered at x = with the same volume (2^)^^ (i.e. ra- 
dius R = (6/7r)i/3£) as that of the cubic region in which 
the function ^(x) in Eq. H66|l increases from —1 to 0. By 
using this expression in Eq. (|51|l . and for q <C i'^, we 
obtain 

A(q)^exp [2 (79)2 In (79)] ^ ^ 

with 7 — {'K/6)^/^i^^ and e{q) = (^qf ■ By studying the 
inverse FT leading from A(q) to -P(F), and therefore to 
W{F), we obtain at large F: 



W{F) ^ 2. (^) 



1/3 



,-4^-3^ 



(69) 



B. Large F prediction for IS. < i/2 



For (5 < 1/2 we can say that the SL falls in the class 
of Sect. IVII A 31 with /3 = 1. It is simple to verify from 
Eq. that ^(x) = — 1 identically in the cube \xi\ < {i— 
2A) with i — 1,2, 3. As written in the previous section, 
in this case P{F) is expected to be rapidly decreasing 
at large F and with finite average quadratic force ap = 
{F^'j as all the higher order moments. 
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More precisely, for 6 -^i 1/2, P(F) is expected to be 
given, to a good approximation, by Eq. H57|) (i.e. Gaus- 
sian with a small kurtosis excess) with a variance ap 
given by Eq. (|29|l where = A^. Instead for S approach- 
ing 1/2, i.e., (1 — 26) <C 1, we expect ap to be given, 
once a suitable isotropic approximation is introduced, by 
Eq. H56|l . A large kurtosis excess Xng for this range of 
6 is also expected, implying large deviations from pure 
Gaussianity. It is particularly interesting to study the 
diverging behavior of ap for S (l/2)~ using Eq. 
to test the validity of this approximated analytical result 
through comparison with measures from numerical sim- 



ulations. In order to apply Eq. H56|l to our SL case we 
need an isotropic approximation as good as possible for 
^(x). First of all it is important to note that (i) for x 
just outside the cube \xi\ < {£ — 2A) with i = 1,2, 3, the 
function ^(x) grows linearly, (ii) ^(x) grows up to the 
surface of the cube \xi\ < 2A ~ £ with i = 1,2,3, (ii) 
outside the cube ^(x) is function with the periodicity of 
the lattice with amplitude at most of order (1 — 25) ^ 1, 
and with zero mean on the period (i.e. on the elemen- 
tary cell). These observations, combined with the same 
argument leading to Eq. Ht)8|) . for S — 1/2, permit one to 
approximate ^(x) simply with 



-1 



e'(^) 



1/3 (x-x') 



2A 



- 1 



for X < X* 



for X > X* 



(70) 



where x* — (6/7r)^/^(^ — 2 A) is chosen so that ^'{x) — 

— 1 in the whole sphere around a; = with the same 
volume 8{£ — 2A)^ as the cube where the exact ^(x) = 

— 1 identically, and {6/'k)^^^£ is analogously the radius 
of the sphere with volume 8^^ [i.e. the volume of the 
cube around x = outside of which ^(x) is everywhere 
small and at most of order (1 — 26)]. Clearly this is a 
rough approximation to ^(x). However we will see that it 
permits one to predict both the logarithmic divergence in 
(1 — 2^) of ap and its order of magnitude. In fact by using 
the function (|7n|l as ^(x) in Eq. ifS^ . with uq = , we 
obtain the following logarithmically diverging behavior 
for 5-> (1/2)- 



of a%: 



4 



-47r 



f- 

V6 



1/3 ln(l - 26) 
25£^ 



(71) 



Moreover, as discussed in the previous section, for values 
of 6 such that (1 — 26) <C 1 the kurtosis excess Xng 
is expected to be large, implying a large deviation from 
Gaussianity of ^'(F). This can be seen by using Eq. H7U|) 
in Eq. (|59|l . which gives for the present case: 



400[ln(l - 26)] 



:{l-25)- 



(72) 



It is important to underline that this result is valid when 
Eq. (|70|l is valid, i.e., when 6 is not too far from the value 
1/2. From Eq. (|72() one can simply see that Xng ~ 0.24 
for 6 = 0.4 and Aatg ~ 1.18 for 5 = 0.44. Moreover, 
beyond this value, Xng diverges monotonously as (1 — 
26)^^ for larger 6. The right interpretation of this result 
is given directly by the statistical meaning of the kurtosis: 
by increasing 6 in this range, W{F) becomes more peaked 
and with a large F tail fatter and fatter than the 3D 
Gaussian distribution H57|l with the same variance a'^. 



C. Large F prediction for A > i!/2 

Since the approximated Chandrasekahr approach to 
the SL problem is more and more precise when 6 in- 
creases, we expect for 6 > 1/2 a better quantitative 
agreement between the analytic results and the numer- 
ical simulations than in the previous two cases. In this 
range of displacements one has ^(0) > —1 meaning that 
on average a particle sees a density of other particles 
larger than zero at a vanishing separation, as each par- 
ticle can be found arbitrarily close at least to another 
particle. This implies that such a SL falls in the class of 
Sect. IVII Ail with a large F tail of W{F) with the same 
scaling as that of the Holtzmark distribution but with a 
different amphtude depending on the value of ^(0). This 
is clear from Eq. H52|l which we rewrite here for conve- 
nience as 



W{F) = [l + m]WH{F) 



(73) 



at large F where WniF) is given by Eq. (|34b|) . As ex- 
plained in Sect. IVII Al as in the homogeneous Poisson 
case, the statistically dominant contribution to the force 
acting on a particle in this cease comes from its nearest 
neighbor. From Eq. H65() it is simple to find that: 



^+2A A 



H 3 



Therefore depending on the choice of A we can have both 
[1 + ^(0)] smaller or larger than one, i.e., with a statistical 
weight for large values of F respectively smaller or larger 
than in the homogeneous Poisson particle distribution, 
depending on whether the probability of finding the near- 
est neighbor at very small distances is smaller or larger 
than in the Poisson case. However as the initial lattice 
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configuration presents negative density-density correla- 
tions at small scales, for most of the choices of A one has 
-1 < ^(0) < and therefore the large F tail of W{F) has 
a smaller amplitude than Wh{F). In general, by taking 
only the largest terms beyond one in Eq. 174(1 . we can 
say that for large A, the large F ratio of W{F)/WHiF) 
approaches unity as [1 — 0(6^"^)] where 0{x) ~ x. 



D. Small F predictions 

For what concerns the small F behavior of W{F) we 
have already seen in the previous section that for all val- 
ues of 6 we have W{F) ~ 47rP(0)-F^ with the pre-factor 
P(0) depending on ^{x). More precisely, for S <^ 1/2 we 
have just seen that P(F) coincides to a good approxima- 
tion with the Gaussian ((57|l . Therefore one has simply 

P(0) = [3/(27rcr|,)] where aj, is given by Eq. For 
higher values of S, when the approximate Chandrasekhar 
approach starts to work, in order to find P{0), one should 
solve the integral H60|l where ^(x) is some appropriate 
isotropic approximation of Eq. ((65|l . Clearly this is a 
task which it is very difficult or impossible to perform 
analytically. However for S > 1/2, i.e., when W{F) is 
power law at large F, one can adopt the following simple 
method to have a rough approximation for the amplitude 
of the small F tail of W{F) and therefore obtain a useful 
approximation of W{F) for all values of F to be used 
to evaluate averages of arbitrary functions of F. One 
assumes the following simple shape for W{F): 

AF^ for F <Fo 
BF-°' for F>Fo, 



W{F) 



(75) 



where respectively, as found above, a = 3 and B = 
27r(7r/6)i/3£-'* for S = 1/2, and a = 5/2 and B = 
27r[l -H^(0)]£-3 [where we have used Eqs. ^ and (|Mall ] 
with ^(0) given by Eq. |(73|) for S > 1/2. In order to find 
A and Fq we impose the following two conditions: (i) 
small F and large F tails take the same value at F — Fq, 
(ii) normalization ofW{F), i.e., dFW{F) = 1. The 
first condition implies 



AF^ 



BFn 



while the second (normalization) condition gives the 
equation: 



A 



B 



-F^ 



= 1 . 



3 " a- 1 
These two equations can be solved to give: 

r A 



BF,-'- 



Fo 



3(a-l) 



3(a+2) 



1/(1-") 



(76) 



In particular it is important to note that for 6 > 1/2 
these formulas imply that 

>o~[i + e(o)]2/3 
A^[i+m]-\ 

where again ^(0) is given by Eq. (|74|l . 



(77) 



E. Comparison with numerical simulations 

To test the above analytical results, we have generated 
numerically several simple cubic SL, with fixed £ (for sim- 
plicity we have chosen £ = I, i.e., no = 1) and different 
values of A in order to study ap and W{F) in a wide 
range going from ^ ^ 1/2 to ^ > 1/2. We expect then to 
see the transition of P(F) from nearly Gaussian for small 
6 to nearly Holtzmark for large values of S when the par- 
ticle distribution approaches the Poisson one. For each 
chosen value of S we have evaluated the PDF W{F) in 
the following way: for each realization of the SL the force 
is evaluated on the "central" particle (i.e. on the particle 
farthest from the boundaries of the system); then W{F) 
is evaluated as a normalized histogram over 10^ realiza- 
tions. The force F on the central particle is computed 
by using the Ewald sum method for lattice sums for the 
cases 6 < 0.3 (i.e. when the SL keeps clear lattice fea- 
tures) in order to make this evaluation faster and precise. 
For larger values of S, on the other hand, F is given by 
the simple sum of the contributions coming from all other 
particles included in the largest sphere centered on the 
central particle. 

In Fig. 2] we present the numerical results for ap vs. 
6 for 6 < 1/2 compared with the theoretical prediction 
for small displacements given by Eqs. H29|l and H3U|) . The 
agreement is excellent up to 6 ~ 0.2. Beyond this value 
ap increases faster than the theoretical prediction for 
small displacements S ^ 1/2, and starts to show the di- 
verging behavior for 5 — > (l/2)~ as predicted by Eqs. H56|) 
and lf7T|l . 

This point is shown better by Fig. |S1 where, in order to 
show the logarithmically diverging behavior of ap when 
S (l/2)~, we have plotted the numerical results for ap 
vs. (1 — 2(5), with a logarithmic scale for the latter, for 
0.4 < 5 < 0.499 (and choosing as above £ = Uq = 1). In- 
deed, if the approximated theoretical prediction Eq. I|71|) 
of a logarithmic divergence is right, this should give a 
straight line. This prediction is verified by the numer- 
ical simulations, albeit with a pre-factor in the loga- 
rithm which is smaller than the theoretical one. This dis- 
crepancy can be explained by the strong approximations 
adopted in Sect.|Vnlto obtain Eq. ||7TIl- **** In Fig.Elwe 
report the comparison for the PDF W{F) between the 
Gaussian theoretical prediction Eq. lf57|l and ap as given 
by the linear approximation (iF^'^p) of Eq. (|29|l for an 
example of the case S <C 1/2 (we have chosen S = 0.05) 
and the numerical results. 

In Fig. [7| on the other hand, we report the numerical 
evaluation of W{F) for the cases S = 0.4 and S — 0.45 vs. 
the Gaussian Wg(F) PDFs with the same variances. As 
theoretically predicted by Eq. (|72|) . already for 5 — 0.4 
it starts to be evident that the actual W{F) has a fatter 
large F tail and a peak lower than that of Wg{F). In 
fact for such values of 5 the kurtosis excess Xng acquires 
a significantly positive value. This discrepancy becomes 
even more clear for S — 0.45 for which Xng ~ 2.1. In fact 
in this case the large F tail starts to develop a power law 
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FIG. 4: Behavior of ap vs. 5 = A/^, for p(u) given by 
Eqs. lHU and ^ with 5 < 1/2. For (5 < 0.2 Eqs. ^ and 
(I3UI I (dashed straight hne), which are vahd only for 5 <^ 1/2, 
applies very well and the agreement with numerical results 
(circles) is excellent up to approximately 5 ~ 0.2. For larger 
values the actual values of ap starts to increase faster than 
this simple linear prediction and the approximation d la Chan- 
drasekhar given in Sect. I VlTl starts to work as shown well by 
the next figure. 
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FIG. 5: Numerical results (circles, and best fit given by the 
dashed line) vs. the approximated theoretical prediction (con- 
tinuous line) Eq. H71|l obtained in Sect. lVIll with £ = no = 1. 
The numerical results, as predicted by Eq. I7H . show a log- 
arithmic divergence of a% in (1 — 2d). However the slope 
is about 20% smaller than the approximated theoretical pre- 
diction. Considering that this theoretical result is obtained 
by a strong approximation (which becomes accurate only for 
larger values of S) , which consists in mimicking the SL with a 
inhomogeneous Poisson particle distribution with radial den- 
sity equal to the average conditional density in the SL, we 
consider this to be a good result. 



FIG. 6: W{F) found numerically for a SL with 5 = 0.05 < 
1/2 and the theoretically predicted Gaussian distribution 
Wg{F) = 4ttF^Pg{F) where Pg(F) is defined by Eq. §7^ 

with ap given by (^IF''^^^ of Eq. The agreement be- 

tween the two curves is very good. 




FIG. 7: Comparison between W{F) found numerically for a 
SL respectively with S = 0.4 and 5 = 0.45 and the Gaussian 
distributions Wg{F) = inF^ Pg{F) with the same variances 
(F^). One observes that in both cases W{F) has its peak 
at smaller F and a more significant large F tail than the 
Gaussian approximation as predicted by Eq. (1721 giving a well 
defined positive kurtosis excess Xng- The closer 5 approaches 
to the "critical" value 1/2, the larger is this deviation. 



feature even though an exponential cut-ofF is still evident. 

The "critical" case S = 1/2 (i.e. where for the first time 
W{F) develops a power law large F tail) is represented 
in Fig. |S| where the numerical W{F) is compared with 
the theoretical prediction for the large and small F tails 
respectively given by Eqs. (|69|) and (|76|l . Despite the 
roughness of the approximation, notably for the small F 
amplitude A, the agreement is very good. 
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FIG. 8: Numerical W{F) from simulations of computer real- 
izations of the SL with 5 = 1/2. The large and small F power 
law approximations are given respectively by Eqs. (1691 and 



FIG. 9: Comparison between numerical W{F) for the case S = 
1 with the exact Holtzmark distribution and the theoretical 
predictions given by Eq. 11731 . with ^(0) given by Eq. (1741 . 
and Eq. ffSl . 



The Holtzmark-like case 6 > 1/2 for W{F) is repre- 
sented in Fig. for the particular value 5 = 1 of the 
shuffling parameter. It is compared both with the ex- 
act Holtzmark distribution obtained in a Poisson particle 
distribution with the same average number density, and 
with the theoretical predictions for the large and small 
F tails given respectively by Eq. H73|) . with ^(0) given by 
Eq. H74|l . and Eq. H76() . On the one hand we see that the 
W{F) approximates quite well the exact Holtzmark one, 
confirming that the shape of W{F) is mainly determined 
by the small separation properties of the particle distribu- 
tion. On the other hand we see also that our theoretical 
approximation shows a good agreement with simulations, 
although the small F prediction is rougher than the large 
F one. This is due to the very simple method we have 
adopted in Sect. I VIII Ul to evaluate the amplitude of this 
tail instead of calculating the more precise but difficult 
Eq. 

Finally, as a further test of our theoretical predictions, 
we have plotted the ratio W{F)/Wh{F) giving a mea- 
sure of the dependence on S of the large F tail of W{F) 
for a wide range of values i5 > 1/2. We have com- 
pared these values with the theoretical prediction given 
by Eq. (O, with ^(0) given by Eq. as functions of S. 
The agreement between numerical simulations and the- 
ory for this quantity is impressive, particularly so given 
the non- monotonous behavior of ^(0). 



IX. DISCUSSION AND CONCLUSIONS 

In this paper we have presented a detailed study of 
the statistical distribution of the total gravitational (or 
Coulomb) force acting on a particle belonging to a ran- 
domly perturbed lattice and due to the sum of the pair 
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FIG. 10: Plot of the W{F)/WHiF) for a wide range of values 
S > 1/2 compared with the theoretical prediction given by 
Eq. with ^(0) given by Eq. 1741 as functions of S. 



gravitational interactions with all the other particles. 

In the first part of the paper we have studied the case 
in which the displacements applied to the lattice particles 
to produce the perturbed lattice are small. In particular 
we have analyzed the linear expansion of the force in the 
displacements. Then we have seen that if only displace- 
ments strictly smaller than half the lattice cell size are 
permitted, this linear expansion can be used to calculate 
to a good approximation the force variance. Otherwise 
this variance goes to infinity, due to the small scale di- 
vergence of the pair interaction, even though the average 
quadratic displacement is kept small. We have seen that 
in the case in which the force variance is finite, it can be 
seen as the sum of two different terms: the former comes 
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from the small random displacements from the lattice 
position of the sources keeping the particle on which the 
force is calculated fixed at its initial lattice position, and 
the latter comes from the displacement of this particle 
from the lattice position keeping at the same time the 
sources at the initial lattice positions. This second term 
can also be seen as the contribution of the uniform nega- 
tive background if it is obtained by summing the contri- 
bution with respect to the original lattice position of the 
particle feeling the force. 

In the subsequent sections we have focused our atten- 
tion on an approximate extension to our case of the Chan- 
drasekhar approach leading to the Holtzmark PDF for 
the homogeneous Poisson particle distribution. In this 
way we have been able to find approximate expressions 
both for the force PDF and its characteristic function 
(and the cumulant generating function) for all the range 
of typical displacements. We have seen that from a qual- 
itative point of view this functional prediction holds for 
the whole range of random displacements, and that the 
agreement becomes quantitatively good for typical dis- 
placements of the order of or larger than half the lattice 
cell size (i.e. when density-density correlations starts to 
be small and the contribution of the Bragg peaks to the 
particle PS is strongly reduced). 

All the above results have then been positively con- 
firmed by a direct comparison with numerical simulation 
of the system in which the SL particle distributions are 
generated with a Monte Carlo like method and the force 
probability distribution is numerically computed. 

We have underlined that, in general, when S starts to 
be of the order of half the lattice cell size, i.e., when the 
minimal permitted distance between particles vanishes, 
the force PDF is dominated by the first NN contribu- 
tion becoming very similar to the Holtzmark distribution 
Wh{F) even though at large distances the SL particle 
distribution is still very different from a homogeneous 
Poisson particle distribution. As has been noted, this is 
due to the small scale divergence of the gravitational (or 
Coulomb) pair interaction between particles. This sug- 
gests that, when the minimal permitted distance between 
particles vanishes, the same behavior for W{F), both at 
small and large F, is expected to be found in all the spa- 
tial particle distributions sharing the same small scale 
correlation properties independently of the large scale 
features. This can be seen clearly in Eq. (|77|l . where 
it is shown that while the small and large F exponent of 
W{F) are universal the amplitudes depend only on ^(0). 

To conclude let us finally return briefly to comment on 
the applications of the results and methods we have just 
found. They can be useful in various different contexts 
mentioned in the introduction, but we will discuss here 
only the primary application which has motivated our 
own study: the comprehension of the dynamics of self- 
gravitating systems studied in cosmology. In this con- 
text large numerical "N-body" simulations of purely self- 
gravitating, essentially point-like |40|] particles are used 
to model the evolution of a self-gravitating fluid. The 



probability distribution of the force on a given particle 
is a useful quantity to understand notably in considering 
(i) the early time dynamics and, more specifically, (ii) 
questions concerning the effects of discreteness in these 
simulations (see 0, [^^l). In a forthcoming work [s^l 
we will report a full analysis of the dynamics of gravita- 
tional evolution of N-body simulations from precisely the 
SL initial conditions analysed here, and the results given 
here will be directly applied in understanding these ques- 
tions. Much can be understood from this study about the 
case of real cosmological N-body simulations, in which 
the initial conditions are lattices subject to small cor- 
related perturbations. In this respect we note that the 
methods developed here, notably the approximate gener- 
alisation of the Chandrasekhar method, can in principal 
be generalised to such distributions. The present study 
of the SL is just a first simpler starting point. Further 
the methods used here can be seen as a first example for 
calculations of other statistical quantities in such distri- 
butions of relevance in understanding the dynamics of 
these self-gravitating systems at larger scales and longer 
times, e.g., the probability distribution of the force on 
the centre of mass of coarse-grained cells, the two-point 
correlation functions of gravitational force as a function 
of separation etc.. 
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APPENDIX A: SMALL q ANALYSIS OF THE 
CUMULANT GENERATING FUNCTION g{q) 

In this appendix we provide some details of the small 
q expansion of the cumulant generating function Q (q) as 
defined in Eqs. ifKHjl and (|5T|l for A > ^/2. As explained 
in Sect. lVlil for such values of A and sufficiently small x 
the correlation function ^(x) can be written as 

ax) = m+Bx^ +oix^), (Ai) 

where ^(0) = -1 for A = £/2 and -1 < ^(0) < -l-oo 
for A > i/2, and B, (3 > 0. Let us suppose that the 
expansion ljAl|l is valid for x < Xq with xq > 0, and 
rewrite the integral in Eq. (|51|l as the following sum of 
two integrals: 
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^(q) 



1 x.fq 
1 sm — 

q \x^ 



+ CXD 



dxx'^[l + £,{x) 



1 ^ ■ ( 1 
1 sm — 



(A2) 



Now in the first integral we can use Eq. IjAip . while in the 
second one, assuming q x^ can expand sin(g/a;^) 
in Taylor series. Since xo > and independent of q, and 
S,{x) vanishes for x +oo, it is now simple to show that 
the second integral is of order q^ at small q. Let us call 
I{q) the first integral, i.e., 



dxx^[l+S,{x)] 



■ { q\ 
— sm — 

q \x'^ J 



By using Eq. (jAl|) we have: 



I{q) = [l+m)]q'^' / dtf [l-t'sinit-^)] +i3g(3+/3)/2 / ^^^2+^ [l - sin^^)] + o 

Jo Jo 



,(3+/3)/2 



, (A3) 



r 



where we have changed variable to i = x/^/q in both 
integrals. In Eq. (|A3|I . the first integral converges in the 
limit q while the second one converges only for < 
(3 < 1, diverges logarithmically for /? = 1, and diverges 



as '^^/^ for /3 > 1. Therefore we can conclude that 
for ^(0) > -1, i.e., A > e/2, in the limit q ^ up to the 
dominant term we have: 



-47rno/(<z) ~ - |4™o[l + C(0)] ^ dtt^ [l - sin(i-2)] | ^3/2 ^ ^ m]CH{q) ■ (A4) 



Instead for ^(0) = —1, i.e., A = 1/2, the coefficient of Note that even for /3 > 1, differently from the case A < 

the first integral of I{q) vanishes, and the second one i'/2 the small g expansion of C/(q) contains a singular part 

[considering also the second integral of Eq. ljA2|l for (3 > even though of order higher than q^ . 
1, which is also of order q^] gives: 

_q(3+/3)/2 for0</3<l 
Q{q)(x { q^\nq for /3 = 1 (A5) 
for P>1. 
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